//Udacity HW 4
//Radix Sorting

#include "reference_calc.cpp"
#include "utils.h"
#include "stdio.h"
/* Red Eye Removal
   ===============
   
   For this assignment we are implementing red eye removal.  This is
   accomplished by first creating a score for every pixel that tells us how
   likely it is to be a red eye pixel.  We have already done this for you - you
   are receiving the scores and need to sort them in ascending order so that we
   know which pixels to alter to remove the red eye.

   Note: ascending order == smallest to largest

   Each score is associated with a position, when you sort the scores, you must
   also move the positions accordingly.

   Implementing Parallel Radix Sort with CUDA
   ==========================================

   The basic idea is to construct a histogram on each pass of how many of each
   "digit" there are.   Then we scan this histogram so that we know where to put
   the output of each digit.  For example, the first 1 must come after all the
   0s so we have to know how many 0s there are to be able to start moving 1s
   into the correct position.

   1) Histogram of the number of occurrences of each digit
   2) Exclusive Prefix Sum of Histogram
   3) Determine relative offset of each digit
        For example [0 0 1 1 0 0 1]
                ->  [0 1 0 1 2 3 2]
   4) Combine the results of steps 2 & 3 to determine the final
      output location for each element and move it there

   LSB Radix sort is an out-of-place sort and you will need to ping-pong values
   between the input and output buffers we have provided.  Make sure the final
   sorted results end up in the output buffer!  Hint: You may need to do a copy
   at the end.

 */

__global__ void histrogram(unsigned int* d_inputValues,unsigned int* d_histroBin,unsigned int* d_LSBout,unsigned int* d_LSBout2,unsigned int mask,int shift_pos,size_t numElems){
    int global_index = threadIdx.x + blockDim.x*blockIdx.x;
    unsigned int bin;
    if(global_index==0){
        d_histroBin[0]=0;
        d_histroBin[1]=0;
    }
    if(global_index < numElems){
        bin = (d_inputValues[global_index] & mask) >> shift_pos;
        atomicAdd(&d_histroBin[!bin],1);
        d_LSBout[global_index]=!bin;
        d_LSBout2[global_index]=bin;
    }
}

__global__ void in_scanblock(unsigned int* d_LSBout,int step){
    unsigned int global_index = threadIdx.x + blockDim.x*blockIdx.x;
    unsigned int ans,step_h;
    extern __shared__ unsigned int s_LSBout[];
    
    s_LSBout[threadIdx.x]=d_LSBout[global_index];
    __syncthreads();
    
    for(int i=0;i<=step;i++){
        step_h = 1<<i;
        if(threadIdx.x + step_h < blockDim.x){
            ans = s_LSBout[threadIdx.x] + s_LSBout[threadIdx.x+step_h];
            __syncthreads();
            s_LSBout[threadIdx.x+step_h] = ans; 
        }
        __syncthreads();
    }
    __syncthreads();
   
    if(threadIdx.x==0){
        for(int j=0;j<blockDim.x;j++){
            d_LSBout[global_index+j]=s_LSBout[j];
        }
    }
}
__global__ void in_scantotal (unsigned int* d_LSBout,unsigned int start_pos){
    d_LSBout[threadIdx.x + start_pos] = d_LSBout[threadIdx.x + start_pos] + d_LSBout[start_pos-1];
}
__global__ void inscan_to_exscan (unsigned int* d_LSBout,size_t numElems){
    unsigned int global_index = threadIdx.x + blockDim.x*blockIdx.x;
    unsigned int temp;
    if(global_index<numElems-1){
        temp=d_LSBout[global_index];
        __syncthreads();
        d_LSBout[global_index+1]=temp;
        d_LSBout[0]=0;
    }

}
__global__ void move_pos(unsigned int* d_inputValues,unsigned int* d_inputPositions,
                         unsigned int* d_LSBout,size_t numElems,unsigned int start_pos,
                         unsigned int* d_outputValues,unsigned int* d_outputPositions){
    unsigned int global_index = threadIdx.x + blockDim.x*blockIdx.x;
    
    if(global_index<numElems-1){
        if(d_LSBout[global_index] < d_LSBout[global_index+1]){
            d_outputValues[d_LSBout[global_index]+start_pos] = d_inputValues[global_index];
            d_outputPositions[d_LSBout[global_index]+start_pos] = d_inputPositions[global_index];
        }
    }
    if(global_index==numElems-1){
        if(d_LSBout[global_index] == d_LSBout[global_index-1]){
            d_outputValues[numElems-1] = d_inputValues[global_index];
            d_outputPositions[numElems-1] = d_inputPositions[global_index];
        }
    }
}
void your_sort(unsigned int* const d_inputVals,
               unsigned int* const d_inputPos,
               unsigned int* const d_outputVals,
               unsigned int* const d_outputPos,
               const size_t numElems)
{
    
    //Copy Input 2-input variables to modifiable variables
    unsigned int *d_inputValues,*d_inputPositions;
    cudaMalloc(&d_inputValues,sizeof(unsigned int)*numElems);
    cudaMalloc(&d_inputPositions,sizeof(unsigned int)*numElems);
    cudaMemcpy(d_inputValues,d_inputVals,sizeof(unsigned int)*numElems,cudaMemcpyDeviceToDevice);
    cudaMemcpy(d_inputPositions,d_inputPos,sizeof(unsigned int)*numElems,cudaMemcpyDeviceToDevice);
    
    unsigned int *d_LSBout,*d_LSBout2,*d_histroBin;
    unsigned int mask=1;
    cudaMalloc(&d_LSBout,sizeof(unsigned int)*numElems);
    cudaMalloc(&d_LSBout2,sizeof(unsigned int)*numElems);
    cudaMalloc(&d_histroBin,sizeof(unsigned int)*2);
    
    unsigned int *h_histroBin = (unsigned int*)malloc(sizeof(unsigned int)*2);
    int blockSize = 848;
    int threadCount = 260;
    
    unsigned int *d_outputValues,*d_outputPositions;
    cudaMalloc(&d_outputValues,sizeof(unsigned int)*numElems);
    cudaMalloc(&d_outputPositions,sizeof(unsigned int)*numElems);
    
for(unsigned int i=0;i<8*sizeof(unsigned int);i++){
    mask=1<<i;
    histrogram<<<1024,1024>>>(d_inputValues,d_histroBin,d_LSBout,d_LSBout2,mask,i,numElems);
    cudaMemcpy(h_histroBin,d_histroBin,sizeof(unsigned int)*2,cudaMemcpyDeviceToHost); 
    /*
    //Histrogram Kernal Debugging
    unsigned int *h_histLSBout = (unsigned int*)malloc(sizeof(unsigned int)*numElems);
    cudaMemcpy(h_histLSBout,d_LSBout,sizeof(unsigned int)*numElems,cudaMemcpyDeviceToHost);
    for(int j=0;j<numElems;j++){
        if(h_histLSBout[j]!=0 && h_histLSBout[j]!=1)
          printf("Error Pos : %d  Value : %u \n",j,h_histLSBout[j]);
    }
    printf("0 Bin : %u   1 Bin :  %u\n",h_histroBin[0],h_histroBin[1]);
    */
    
    /*
    //Exclusive Scan Kernal Debugging-1
    
    unsigned int *h_LSBout=(unsigned int*)malloc(sizeof(unsigned int)*numElems);
    unsigned int *h_serial_exscan=(unsigned int*)malloc(sizeof(unsigned int)*numElems);
    unsigned int *h_kernal_inscan=(unsigned int*)malloc(sizeof(unsigned int)*numElems);
    cudaMemcpy(h_LSBout,d_LSBout,sizeof(unsigned int)*numElems,cudaMemcpyDeviceToHost);
    
    h_serial_exscan[0]=0;
    for(int j=1;j<numElems;j++){
        h_serial_exscan[j] = h_serial_exscan[j-1]+h_LSBout[j-1];
    }
    */
    //Total Element : 220480  We split them into 530 Block , 416 Threads per blocks
    
    in_scanblock<<<blockSize,threadCount,sizeof(unsigned int)*threadCount>>>(d_LSBout,ceil(log(threadCount)/log(2)));
    for(int i=threadCount;i<blockSize*threadCount;i=i+threadCount){
        in_scantotal<<<1,threadCount>>>(d_LSBout,i);
    }
    inscan_to_exscan<<<1024,1024>>>(d_LSBout,numElems);
    
    in_scanblock<<<blockSize,threadCount,sizeof(unsigned int)*threadCount>>>(d_LSBout2,ceil(log(threadCount)/log(2)));
    for(int j=threadCount;j<blockSize*threadCount;j=j+threadCount){
        in_scantotal<<<1,threadCount>>>(d_LSBout2,j);
    }
    inscan_to_exscan<<<1024,1024>>>(d_LSBout2,numElems);
    /*
    //Inclusive Scan Kernal Debugging-2
    cudaMemcpy(h_kernal_inscan,d_LSBout,sizeof(unsigned int)*numElems,cudaMemcpyDeviceToHost);
    for(int j=0;j<numElems;j++){
        if(h_kernal_inscan[j]!=h_serial_exscan[j]){
          printf("Error Pos : %d Value : %u Serial: %u  Kernal: %u \n",j,h_LSBout[j],h_serial_exscan[j],h_kernal_inscan[j]);
          break;
        }
    }
    */
    //Move Values&Positions to their pos
    move_pos<<<1024,1024>>>(d_inputValues,d_inputPositions,d_LSBout,numElems,0,d_outputValues,d_outputPositions);
    move_pos<<<1024,1024>>>(d_inputValues,d_inputPositions,d_LSBout2,numElems,h_histroBin[1],d_outputValues,d_outputPositions);
    /*
    //Move Position Debugging
    unsigned int *h_outputValues = (unsigned int*)malloc(sizeof(unsigned int)*numElems);
    cudaMemcpy(h_outputValues,d_outputValues,sizeof(unsigned int)*numElems,cudaMemcpyDeviceToHost);
    for(int j=0;j<numElems;j++){
        if((h_outputValues[j]&1)!=1 && (h_outputValues[j]&1)!=0){
            printf("Values %u , Output %u \n",h_outputValues[j],h_outputValues[j]&1);
            break;
        }
    }
    */
    cudaMemcpy(d_inputValues,d_outputValues,sizeof(unsigned int)*numElems,cudaMemcpyDeviceToDevice);
    cudaMemcpy(d_inputPositions,d_outputPositions,sizeof(unsigned int)*numElems,cudaMemcpyDeviceToDevice);
    
    unsigned int *h_outputval=(unsigned int*)malloc(sizeof(unsigned int)*numElems);
    unsigned int *h_outputpositions=(unsigned int*)malloc(sizeof(unsigned int)*numElems);
    cudaMemcpy(h_outputval,d_outputValues,sizeof(unsigned int)*numElems,cudaMemcpyDeviceToHost);
    cudaMemcpy(h_outputpositions,d_outputPositions,sizeof(unsigned int)*numElems,cudaMemcpyDeviceToHost);
    for(int j=0;j<20;j++){
        printf("Position: %u Value:%u \n",h_outputpositions[j],h_outputval[j]);
    }
    printf("\n**********************\n");
}
    cudaMemcpy(d_outputVals,d_outputValues,sizeof(unsigned int)*numElems,cudaMemcpyDeviceToDevice);
    cudaMemcpy(d_outputPos,d_outputPositions,sizeof(unsigned int)*numElems,cudaMemcpyDeviceToDevice);
    
    
    unsigned int *h_inputval=(unsigned int*)malloc(sizeof(unsigned int)*numElems);
    unsigned int *h_positions=(unsigned int*)malloc(sizeof(unsigned int)*numElems);
    cudaMemcpy(h_inputval,d_inputVals,sizeof(unsigned int)*numElems,cudaMemcpyDeviceToHost);
    cudaMemcpy(h_positions,d_inputPos,sizeof(unsigned int)*numElems,cudaMemcpyDeviceToHost);
    
    unsigned int *h_outputval=(unsigned int*)malloc(sizeof(unsigned int)*numElems);
    unsigned int *h_outputpositions=(unsigned int*)malloc(sizeof(unsigned int)*numElems);
    cudaMemcpy(h_outputval,d_outputValues,sizeof(unsigned int)*numElems,cudaMemcpyDeviceToHost);
    cudaMemcpy(h_outputpositions,d_outputPositions,sizeof(unsigned int)*numElems,cudaMemcpyDeviceToHost);
    unsigned int max=h_inputval[0];
    unsigned int min=h_inputval[0];
    for(int i=0;i<numElems;i++){
        if(h_inputval[i] > max){
            max=h_inputval[i];
        }
        if(h_inputval[i]<min){
            min=h_inputval[i];
        }
    }
    printf("Serial : %u  Kernal : %u \n",min,h_outputval[0]);
    printf("Serial : %u  Kernal : %u %u %u %u\n",max,h_outputval[numElems-1],h_outputval[numElems-2],h_outputval[numElems-3],h_outputval[numElems-4]);
    cudaFree(d_inputValues);
    cudaFree(d_inputPositions);
    cudaFree(d_LSBout);
    cudaFree(d_histroBin);
    cudaFree(d_outputValues);
    cudaFree(d_outputPositions);
    //free(h_histroBin);
    //free(h_histLSBout);
    //free(h_LSBout);
    //free(h_serial_exscan);
    //free(h_kernal_inscan);
}